home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / demos / velpic / velpic.sci < prev   
Text File  |  1999-09-16  |  33KB  |  1,116 lines
  1. //[vel,regionlist,linelist,seedlist,velolist]=velpic(nz,nx,sext)
  2. //[vel,regionlist,linelist,seedlist]=velpic(nz,nx)
  3. //Macro which interactively defines regions of a matrix.
  4. //Closed surfaces inside the frame are not allowed.
  5. // nz         :Number of columns in the matrix
  6. // nx         :Number of rows in the matrix
  7. // vel        :velocity matrix
  8. // sext       :[smin,smax] gives the velocities min and max
  9. // regionlist :list whose elements are matrices of dimension
  10. //            :(2xN) containing the indices of the different regions
  11. // linelist   :list whose elements are matrices of dimension (2xM)
  12. //            :containing the indices of the different lines
  13. // seedlist   :(2xK) vector containing positions of seeds
  14. // velolist   :(1xK) vector containing velocities associated to seeds
  15. //
  16. //!
  17. // author: C. Bunks     date: 12-NOV-90
  18. [lhs,rhs]=argn(0)
  19. //turn off the scilab function 'more'
  20.  
  21.    nrnc=lines();
  22.    lines(0);
  23. //
  24.    qt='off'
  25. //
  26.    select rhs
  27.    case 2 then
  28.       smin=1500;smax=6000;
  29.       rnd=100
  30.    case 3 then
  31.       smin=sext(1);smax=sext(2)
  32.       rnd=10**int(-1+log(abs(smax-smin))/log(10))
  33.  
  34.    else error('velpic(nz,nx [,sext])')
  35.    end
  36. //make a suitable frame
  37.  
  38.    bnames=list('stop','help','undo','grill','quit');
  39.    gopt='off';
  40.    [buttons,slides]=makeframe(nz,nx,bnames);
  41.    slide1=slides(1,:);
  42.  
  43. //seperate matrix into regions defined by lines
  44.  
  45.    linelist=makehorizons(nz,nx,bnames,buttons);
  46.    if qt='on' then return,end
  47. //determine indices of matrix in each region
  48.  
  49.    write(%io(2),'Searching Regions'),
  50.    regionlist=list();
  51.    indexlist=[ones(1,nz).*.(1:nx);(1:nz).*.ones(1,nx)];
  52.  
  53.    region=0;
  54.    while size(indexlist)<>[0 0],
  55.       region=region+1;
  56.       seed=indexlist(:,1);
  57.       [inds,indexlist]=id_rgn(indexlist,linelist,seed);
  58.       regionlist(region)=inds;
  59.       write(%io(2),' Region Identified:  '+string(region)),
  60.    end,
  61.  
  62. //sow a seed in each region
  63.  
  64.    [seedlist]=sow(nz,nx,slide1,bnames,buttons,regionlist,linelist);
  65.    if qt='on' then return,end
  66.    [ssl,isl]=sort(seedlist(3,:));//sort seeds and velocities to regions
  67.    [sr,sc]=size(seedlist);
  68.    velolist=seedlist(4,isl(sc:-1:1));
  69.    rlist=seedlist(3,isl(sc:-1:1));
  70.    seedlist=seedlist(1:3,isl(sc:-1:1));
  71.  
  72. //Load velocity matrix with velocities by region.
  73.  
  74.    i2=0;
  75.    vel=0*ones(1,nz*nx);
  76.    nor=size(regionlist);
  77.    nov=maxi(size(velolist));
  78.    write(%io(2),'Building velocity matrix'),
  79.    for k=1:nor,
  80.       i1=i2+1;i2=i1-1;
  81.       for j=1:nov,
  82.          if rlist(j)=k then, i2=i2+1; end,
  83.       end,
  84.       rk=regionlist(k);
  85.       rv=(rk(1,:)-ones(rk(1,:)))*nz+rk(2,:);
  86.       if i1=i2 then,//constant velocity region
  87.          vel(rv)=velolist(i1:i2)*ones(rv);
  88.       else,//linearly varying velocity region
  89.          vel(rv)=velcalc(rk,seedlist(:,i1:i2),velolist(i1:i2));
  90.       end,
  91.    end,
  92.    vel=matrix(vel,nz,nx);
  93.    vel=vel(nz:-1:1,:);
  94.  
  95. //re-establish the line counter
  96.  
  97.    lines(nrnc(2));
  98.  
  99. //end
  100.  
  101. //[]=helpme(msn)
  102. //[]=helpme(msn)
  103. //macro containing help messages
  104. // msn :message number
  105. //
  106. //!
  107. // author: C. Bunks     date: 12-NOV-90
  108.  
  109. hm1=[' ';
  110. 'Begin drawing a line by clicking the left';
  111. 'mouse button across a previously drawn line';
  112. 'segment or across a frame line segment. In';
  113. 'the second case the first mouse click must';
  114. 'be in the margin.';
  115. ' ';
  116. 'STOP:';
  117. '  Stop drawing lines and start picking velocities';
  118. ' ';
  119. 'HELP:';
  120. '  This help message';
  121. ' ';
  122. 'UNDO:';
  123. '  Undo a previously drawn line';
  124. ' ';
  125. 'GRILL:';
  126. '  Toggle grill on/off';
  127. ' '];
  128.  
  129. hm2=[' ';
  130. 'Terminate drawing a line by clicking the left';
  131. 'mouse button across a previously drawn line';
  132. 'segment or across a frame line segment.';
  133. ' ';
  134. 'STOP:';
  135. '  Completely undo the current line and begin';
  136. '  drawing a new line';
  137. ' ';
  138. 'HELP:';
  139. '  This help messaage';
  140. ' ';
  141. 'UNDO:';
  142. '  Undo a line segment';
  143. ' ';
  144. 'GRILL:';
  145. '  Toggle grill on/off';
  146. ' '];
  147.  
  148. hm3=[' ';
  149. 'To sow a seed, choose a velocity by clicking';
  150. 'the left mouse button in the velocity slide.';
  151. ' ';
  152. 'The velocity value chosen on the slide';
  153. 'is echoed to the terminal.';
  154. ' ';
  155. 'To change the velocity value click the left';
  156. 'mouse button at another location on the slide.';
  157. ' ';
  158. 'When a suitable velocity value has been obtained';
  159. 'click the left mouse button somewhere in the frame';
  160. 'to plant the velocity seed';
  161. ' ';
  162. 'Multiple seeds in the same region yields linear';
  163. 'velocity gradients.  The linear gradients are';
  164. 'established between consecutive seed values';
  165. '(following the order of placement)';
  166. ' ';
  167. 'STOP:';
  168. '  Stop planting seeds and return velocity matrix';
  169. ' ';
  170. 'HELP:';
  171. '  This help message';
  172. ' ';
  173. 'UNDO:';
  174. '  Undo a seed placement';
  175. ' ';
  176. 'GRILL:';
  177. '  Toggle grill on/off';
  178. ' '];
  179.  
  180. hmt=list(hm1,hm2,hm3);
  181. x_message(hmt(msn));
  182. //end
  183.  
  184. //[seedlist]=sow(nz,nx,slide,bnames,buttons,regionlist,linelist);
  185. //[seedlist]=sow(nz,nx,slide,bnames,buttons,regionlist,linelist);
  186. //Interactively identify the various regions by clicking
  187. //the mouse in each region
  188. // nz         :number of rows of matrix
  189. // nx         :number of columns of matrix
  190. // nol        :number of lines (not including the boundary)
  191. // slide      :slide information
  192. // bnames     :button names
  193. // buttons    :button positions
  194. // regionlist :list of regions
  195. // seedlist   :(4xN) matrix where each column contains (from top to
  196. //            :bottom) the x and y coordinates of the seed, the region
  197. //            :that the seed is in, and the velocity value associated
  198. //            :to the seed
  199. //
  200. //!
  201. // author: C. Bunks     date: 12-NOV-90
  202.  
  203.    sl1=slide(1);sl2=slide(2);sl3=slide(3);
  204.    sl4=slide(4);sl5=slide(5);sl6=slide(6);
  205.    seedlist=[];
  206.  
  207. //slide bar constants  
  208.  
  209.    vbar=[0 0 0 0 0;0 0 0 0 0];
  210.    if nx<=nz then,
  211.       sw=sl2-sl1;
  212.       b1=sl1+.05*sw;b2=sl2-.05*sw;
  213.       b3=.05*sw;b4=.05*sw;
  214.    else,
  215.       sw=sl4-sl3;
  216.       b1=.05*sw;b2=.05*sw
  217.       b3=sl3+.05*sw;b4=sl4-.05*sw;
  218.    end,
  219.  
  220. //get velocity and place seed
  221.  
  222.    nol=size(linelist);
  223.    nor=size(regionlist);
  224.  
  225.    sflag='on';
  226.    while sflag='on',
  227.       write(%io(2),'Choose Velocity'),
  228.       veloflag='on';
  229.       vel=-1;
  230.       while veloflag='on',//loop until desired vel is obtained
  231.           [i_i,v1,v2]=xclick();v=[v1;v2];
  232.  
  233. //check for a button
  234.  
  235.    [br,bc]=size(buttons);
  236.    hflag='on';
  237.    while hflag='on',
  238.       hflag='off';
  239.       for bi=1:br,
  240.          if buttons(bi,1)<v1 then, if v1<buttons(bi,2) then,
  241.          if buttons(bi,3)<v2 then, if v2<buttons(bi,4) then,
  242.             select bnames(bi)
  243.             case 'help' then,
  244.                helpme(3);
  245.                [i_i,v1,v2]=xclick();v=[v1;v2];
  246.                hflag='on';
  247.             case 'stop' then,//stop sowing (all regions >= one seed)
  248.                [sr,sc]=size(seedlist);
  249.                rflag='off';
  250.                for rk=1:nol,
  251.                   srflag='off';
  252.                   for sk=1:sc,
  253.                      if seedlist(3,sk)=rk then, srflag='on'; end,
  254.                   end,
  255.                   if srflag='off' then, rflag='on'; end,
  256.                end,
  257.                if rflag='on' then,
  258.                   write(%io(2),'Not every region has a seed'),
  259.                   write(%io(2),'Choose Velocity'),
  260.                   [i,v1,v2]=xclick();v=[v1;v2];
  261.                   hflag='on';
  262.                else,
  263.                   sflag='off';
  264.                   veloflag='off';
  265.                end,
  266.             case 'grill' then,
  267.                if gopt='on' then, 
  268.                   gopt='off';
  269.                else,
  270.                   gopt='on';
  271.                end,
  272.                makegrill(nx,nz,gopt);
  273.                [i,v1,v2]=xclick();v=[v1;v2];
  274.                hflag='on';
  275.             case 'undo' then,//undo line segment
  276.                [sr,sc]=size(seedlist);
  277.                toff=.05*maxi([nx,nz])/3;
  278.                if sc>0 then,
  279.                   s1=seedlist(1,sc);
  280.                   s2=seedlist(2,sc);
  281.                   text=string(seedlist(4,sc));
  282.                   dess(36)=15;
  283.                   xset("alufunction",6);
  284.                     plot2d(s1',s2',[3,0],"000");
  285.                     xstring(s1+toff,s2+toff,text,0,0);
  286.                   xset("alufunction",3);
  287.                   seedlist=seedlist(:,1:sc-1);
  288.                else,
  289.                   write(%io(2),'Nothing to undo'),
  290.                   write(%io(2),'Choose Velocity'),
  291.                end,
  292.                [i,v1,v2]=xclick();v=[v1;v2];
  293.                hflag='on';
  294.              case 'quit' then,
  295.                qt=resume('on');
  296.             end,
  297.          end,end,
  298.          end,end,
  299.       end,
  300.    end,
  301.  
  302.       if sflag='on' then,//stop hasn't been signalled
  303. //get velocity
  304.  
  305.          if sl1<=v1 then, if v1<=sl2 then,
  306.          if sl3<=v2 then, if v2<=sl4 then,
  307. //calculate velocity rounded to nearest rnd
  308. //plot bar indicator and give velocity value
  309.             xset("alufunction",6);           
  310.             plot2d(vbar(1,:)',vbar(2,:)',[-1],"000");
  311.             xset("alufunction",3);           
  312.             if nx<=nz then,
  313.                vel=rnd*round((sl5+(sl6-sl5)*(v2-sl3)/(sl4-sl3))/rnd);
  314.                vbar=[b1 b2 b2 b1 b1;v2-b3 v2-b3 v2+b4 v2+b4 v2-b3];
  315.                plot2d(vbar(1,:)',vbar(2,:)',[-1],"000"),
  316.             else,
  317.                vel=rnd*round((sl5+(sl6-sl5)*(v1-sl1)/(sl2-sl1))/rnd);
  318.                vbar=[v1-b1 v1+b2 v1+b2 v1-b1 v1-b1;b3 b3 b4 b4 b3];
  319.                plot2d(vbar(1,:)',vbar(2,:)',[-1],"000"),
  320.             end,
  321.             write(%io(2),vel,'(f10)'),
  322.          end,end,
  323.          end,end,
  324.          if vel<>-1 then,//check that a velocity was chosen
  325. //plant seed
  326.             if 1<=v1 then, if v1<=nx then,
  327.             if 1<=v2 then, if v2<=nz then,
  328.  
  329. //find the region that v is in
  330.  
  331.                rflag='on';
  332.                nr=0;
  333.                while rflag='on',//look for a region that contains seed
  334.                   nr=nr+1;
  335.                   rk=regionlist(nr);
  336.                   rflag='off';
  337.                   for nl=1:nol,//for this region check all lines
  338.                      [testflag,bav]=testpt(v,rk(:,1),linelist(nl));
  339.                      if testflag='on' then, rflag='on'; end,
  340.                   end,
  341.                   if rflag='off' then,//this region if no intersections
  342.                      sregion=nr;   
  343.                   end,
  344.                end,
  345.  
  346. //plot seed and add to seedlist
  347.                mrgn=.05*maxi([nx,nz])/3;
  348.                plot2d(v1',v2',[3,0],"000");
  349.                xstring(v1+mrgn,v2+mrgn,string(vel),0,0);
  350.                seedlist=[seedlist,[v;sregion;vel]];
  351.                veloflag='off';
  352.             end,end,
  353.             end,end,
  354.          end,
  355.       end,
  356.       end,
  357.    end,
  358.  
  359. //end
  360.  
  361. //[linelist]=makehorizons(nz,nx,bnames,buttons)
  362. //[linelist]=makehorizons(nz,nx,bnames,buttons)
  363. //macro which creates a frame with control buttons and
  364. //interactively draws lines (i.e., horizons)
  365. // nz       :Number of rows in frame
  366. // nx       :Number of columns in frame
  367. // bnames   :Button names
  368. // buttons  :Button locations
  369. // linelist :list whose elements are (2xN) matrices of line indices
  370. //
  371. //!
  372. // author: C. Bunks     date: 12-NOV-90
  373.  
  374. //define outer perimeter as a line
  375.  
  376.    prt=.001;
  377.    linelist=list([1-prt 1-prt  nx+prt nx+prt 1-prt;...
  378.                   1-prt nz+prt nz+prt 1-prt  1-prt]);
  379.    nol=size(linelist);
  380.  
  381. //start line drawing
  382.  
  383.    layer='true';
  384.    yec=[];
  385.    while layer='true',
  386.       nol=nol+1;
  387.  
  388. //Define layer by drawing a line
  389.  
  390.       write(%io(2),'Draw a Line'),
  391.       [line,linelist,yec]=drawline(nz,nx,linelist,yec,bnames,buttons);
  392.       if size(line)=[0,0] then,//if returned line is empty then stop
  393.          layer='false'; 
  394.       else if line=0 then,
  395.          nol=nol-1;
  396.       else,
  397.          nol=size(linelist);
  398.          nol=nol+1;
  399.          linelist(nol)=line;
  400.       end,
  401.       end,
  402.    end,
  403.  
  404. //end
  405.  
  406. //[buttons,slides]=makeframe(nz,nx,btextlist);
  407. //[buttons,slides]=makeframe(nz,nx,btextlist);
  408. //macro which makes a frame and buttons
  409. // nz        :Number of rows in frame
  410. // nx        :Number of columns in frame
  411. // btextlist :text for buttons
  412. // buttons   :button information
  413. // slides    :slide information
  414. //
  415. //!
  416. // author: C. Bunks     date: 12-NOV-90
  417.  
  418. //setup of frame
  419.  
  420. //draw work box.
  421. //note that the work box is placed in [1,nx] x [1,nz]
  422. //and that the remaining elements are for centering properly
  423. //the box (i.e., mrgn and the 1).
  424.  
  425.    max=maxi([nz,nx]);
  426.    mrgn=.05*(max-1);
  427.    xmin=1-mrgn-max+nx;xmax=max+4*mrgn;
  428.    ymin=1-mrgn-max+nz;ymax=max+4*mrgn;
  429.    if nx>nz+4*mrgn then,
  430.       xmin=1-mrgn-max+nx;xmax=max+mrgn;
  431.       ymin=1-mrgn-max+nz;ymax=max+mrgn;
  432.    end,
  433.    if nx<nz-4*mrgn
  434.       xmin=1-mrgn-max+nx;xmax=max+4*mrgn;
  435.       ymin=1-mrgn-max+nz;ymax=max+mrgn;
  436.    end,
  437.    rect=[xmin,ymin,xmax,ymax];
  438.    plot2d(0,0,[-1],"012",' ',rect);
  439.    plot2d([1;1;nx;nx;1],[nz;1;1;nz;nz],[-1],"000"),
  440. //make buttons and slides
  441.    dess4=10
  442.    dess6=10
  443.    dx=(xmax-xmin)/dess4
  444.    dy=(ymax-ymin)/dess6
  445.  
  446. //make button bank
  447.  
  448.    bs=size(btextlist);//number of buttons
  449.    bsi=1/(2*bs);
  450.    buttons=[];
  451.    if nz=>nx then,//in the right side margin
  452.       bc=int(bs/(4+bsi))+1;//number of columns
  453.       br=int(bs/(bc+bsi))+1;//number of rows
  454.       bm=2*mrgn/(bc+1);
  455.       bentry=0;
  456.       for bj=1:bc,for bi=1:br,
  457.          xbmin=nx+2*(bj-1)*bm+bm;xbmax=xbmin+2*bm;
  458.          ybmax=nz-2*(bi-1)*bm;ybmin=ybmax-2*bm;
  459.          bentry=bentry+1;
  460.          if bentry<=bs then,
  461.             text=btextlist(bentry);
  462.          else,
  463.             text='    ';
  464.          end,
  465.          makebutton(xbmin,xbmax,ybmin,ybmax,dx,dy,text);
  466.          buttons=[buttons;xbmin,xbmax,ybmin,ybmax];
  467.       end,end,
  468.    else,//in the top margin
  469.       br=int(bs/(4+bsi))+1;//number of rows
  470.       bc=int(bs/(br+bsi))+1;//number of columns
  471.       bm=2*mrgn/(br+1);
  472.       bentry=0;
  473.       for bi=1:br,for bj=1:bc,
  474.          xbmin=1+2*(bj-1)*bm;xbmax=xbmin+2*bm;
  475.          ybmin=nz+2*(bi-1)*bm+bm;ybmax=ybmin+2*bm;
  476.          bentry=bentry+1;
  477.          if bentry<=bs then,
  478.             text=btextlist(bentry);
  479.          else,
  480.             text='    ';
  481.          end,
  482.          makebutton(xbmin,xbmax,ybmin,ybmax,dx,dy,text);
  483.          buttons=[buttons;xbmin,xbmax,ybmin,ybmax];
  484.       end,end,
  485.    end,
  486.  
  487. //velocity slide
  488.  
  489.    if nz=>nx then,
  490.       xbmin=nx+mrgn;xbmax=xbmin+2*mrgn;
  491.       ybmin=1;ybmax=(1+nz)/2;
  492.    else,
  493.       xbmin=(1+nx)/2;xbmax=nx;
  494.       ybmin=nz+mrgn;ybmax=ybmin+2*mrgn;
  495.    end,
  496.    text='VEL';
  497.    dx=(xmax-xmin)/dess4;
  498.    dy=(ymax-ymin)/dess6;
  499.    makeslide(xbmin,xbmax,ybmin,ybmax,dx,dy,text,smin,smax);
  500.    slides=[xbmin,xbmax,ybmin,ybmax,smin,smax];
  501.  
  502. //end
  503.  
  504. //[]=makebutton(xbmin,xbmax,ybmin,ybmax,dx,dy,text)
  505. //[]=makebutton(xbmin,xbmax,ybmin,ybmax,dx,dy,text)
  506. //make to make a button
  507. // xbmin :min x coordinate of button
  508. // xbmax :max x coordinate of button
  509. // ybmin :min y coordinate of button
  510. // ybmax :max y coordinate of button
  511. // dx    :ratio of plot length to frame length: (xmax-xmin)/dess(4)
  512. // dy    :ratio of plot height to frame height: (ymin-ymax)/dess(6)
  513. // text  :text in button
  514. //
  515. //!
  516. // author: C. Bunks     date: 12-NOV-90
  517. //make button box
  518.    xstringb(xbmin,ybmin,text,xbmax-xbmin,ybmax-ybmin);
  519.    xrect(xbmin,ybmax,xbmax-xbmin,ybmax-ybmin);
  520. //end
  521.  
  522. //[]=makeslide(xbmin,xbmax,ybmin,ybmax,dx,dy,text,smin,smax)
  523. //[]=makeslide(xbmin,xbmax,ybmin,ybmax,dx,dy,text)
  524. //make to make a button
  525. // xbmin :min x coordinate of button
  526. // xbmax :max x coordinate of button
  527. // ybmin :min y coordinate of button
  528. // ybmax :max y coordinate of button
  529. // dx    :ratio of plot length to frame length: (xmax-xmin)/dess(4)
  530. // dy    :ratio of plot height to frame height: (ymin-ymax)/dess(6)
  531. // text  :text in button
  532. // smin  :min value of slide
  533. // smax  :max value of slide
  534. //
  535. //!
  536. // author: C. Bunks     date: 12-NOV-90
  537.  
  538. //NOTE: The constant wcf is a 'wierd correction factor necessary
  539. //for the correct positioning of the text in the button.  I don't
  540. //understand why it is needed...
  541.  
  542.    wcf=1.3;
  543.  
  544. //make slide box
  545.  
  546.    xvec=[xbmin xbmin xbmax xbmax xbmin];//x button frame vector
  547.    yvec=[ybmax ybmin ybmin ybmax ybmax];//y button frame vector
  548.    delx=xbmax-xbmin;dely=ybmax-ybmin;
  549.  
  550. //label slide
  551.  
  552.    dess(27)=1;//soft caracters
  553.    sl=length(text);
  554.    if delx<dely then,
  555.       cw=(xbmax-xbmin)/(dx*sl);//char.width=.9*(box width)/(text length)
  556.    else,
  557.       cw=(ybmax-ybmin)/(dx*sl);
  558.    end,
  559.    ch=cw/.82;
  560. //   setcsize(cw,ch);
  561.    textx=(xbmin+xbmax-dx*sl*cw/wcf)/2;//text x position
  562.    texty=(ybmin+ybmax-dy*ch/wcf)/2;//text y position
  563.    dess(51)=2;//turn labelling on
  564.    dess(50)=textx;dess(52)=texty;//positions of text
  565.    plot2d(xvec',yvec',[-1],"000");
  566.    xstring(textx,texty,text,0,0),
  567.  
  568. //label min value of slide
  569.  
  570.    text=string(smin);
  571.    sl=length(text);
  572.    if delx<dely then,
  573.       cw=(xbmax-xbmin)/(dx*sl);
  574.       ch=cw/.82;
  575.       textx=(xbmin+xbmax-dx*sl*cw/wcf)/2;//text x position
  576.       texty=ybmin+dy*ch/wcf/2;//text y position
  577.    else,
  578.       cw=(ybmax-ybmin)/(dx*sl);
  579.       ch=cw/.82;
  580.       textx=xbmin+dx*cw/wcf;//text x position
  581.       texty=(ybmin+ybmax-dy*cw/wcf)/2;//text y position
  582.    end,
  583. //   setcsize(cw,ch);
  584.    dess(51)=2;//turn labelling on
  585.    dess(50)=textx;dess(52)=texty;//positions of text
  586.    xstring(textx,texty,text,0,0),
  587.  
  588. //label max value of slide
  589.  
  590.    text=string(smax);
  591.    sl=length(text);
  592.    if delx<dely then,
  593.       cw=(xbmax-xbmin)/(dx*sl);
  594.       ch=cw/.82;
  595.       textx=(xbmin+xbmax-dx*sl*cw/wcf)/2;//text x position
  596.       texty=ybmax-3*dy*ch/wcf/2;//text y position
  597.    else,
  598.       cw=(ybmax-ybmin)/(dx*sl);
  599.       ch=cw/.82;
  600.       textx=xbmax-dx*(1+sl)*cw/wcf;//text x position
  601.       texty=(ybmin+ybmax-dy*cw/wcf)/2;//text y position
  602.    end,
  603. //   setcsize(cw,ch);
  604.    dess(51)=2;//turn labelling on
  605.    dess(50)=textx;dess(52)=texty;//positions of text
  606.    xstring(textx,texty,text,0,0),
  607.  
  608. //end
  609.  
  610. //[line,linelist,yec]=drawline(nz,nx,linelist,yec,bnames,buttons)
  611. //[line,linelist,yec]=drawline(nz,nx,linelist,yec,bnames,buttons)
  612. //interactively draw a line
  613. // nz       :number of columns of matrix
  614. // nx       :number of rows of matrix
  615. // linelist :list containing the matrices of indices for 
  616. //          :each line
  617. // yec      :extension list
  618. // bnames   :Names of buttons
  619. // buttons  :Button locations
  620. // line     :new line
  621. //
  622. //!
  623. // author: C. Bunks     date: 12-NOV-90
  624.  
  625. //to begin new line the first two clicks of the mouse must
  626. //intersect an old line or the frame
  627.  
  628.    line=[];
  629.    [i,x1,x2]=xclick();
  630.    testflag='off';
  631.    flag='on';
  632.  
  633. //check for a button
  634.  
  635.    [br,bc]=size(buttons);
  636.    hflag='on';
  637.    while hflag='on',//while the mouse has been clicked in a button
  638.       hflag='off';
  639.       for bi=1:br,//check all the buttons
  640.          if buttons(bi,1)<x1 then, if x1<buttons(bi,2) then,
  641.          if buttons(bi,3)<x2 then, if x2<buttons(bi,4) then,
  642.             select bnames(bi)
  643.             case 'help' then,
  644.                helpme(1);
  645.                [i_i,x1,x2]=xclick();
  646.                hflag='on';
  647.             case 'stop' then,
  648.                testflag='on';
  649.                flag='off';
  650.                yecl=0;
  651.             case 'grill' then,
  652.                if gopt='on' then, 
  653.                   gopt='off';
  654.                else,
  655.                   gopt='on';
  656.                end,
  657.                makegrill(nx,nz,gopt);
  658.                [i_i,x1,x2]=xclick();
  659.                hflag='on';
  660.             case 'undo' then,//undo last line
  661.                ls=size(linelist);
  662.                if ls>1 then,//remove a line
  663.                   [yer,yecc]=size(yec);
  664.                   yf=yec(1,yecc);
  665.                   yl=yec(2,yecc);
  666.                   lk=linelist(ls);
  667.                   [lkr,lkc]=size(lk);
  668.                   xset("alufunction",6);
  669.                   plot2d(lk(1,yf+1:lkc-yl)',lk(2,yf+1:lkc-yl)',[-1],"000"),
  670.                   xset("alufunction",3);
  671.                   ltemp=list();
  672.                   for k=1:ls-1,
  673.                      ltemp(k)=linelist(k);
  674.                   end,
  675.                   linelist=ltemp;
  676.                   yec=yec(:,1:yecc-1);
  677.                else,
  678.                   write(%io(2),'Nothing to undo'),
  679.                end,
  680.                [i_i,x1,x2]=xclick();
  681.                hflag='on';
  682.             case 'quit' then,
  683.                qt=resume('on')
  684.             end,
  685.          end,end,
  686.          end,end,
  687.       end,
  688.    end,
  689.  
  690. //start drawing line
  691.  
  692.    yext=[];
  693.    while testflag='off',   
  694.       plot2d(x1',x2',[3,1],"000");//make a start circle
  695.       [i_i,y1,y2]=xclick();
  696.  
  697. //check for a button
  698.  
  699.    [br,bc]=size(buttons);
  700.    hflag='on';
  701.    while hflag='on'
  702.       hflag='off';
  703.       for bi=1:br,
  704.          if buttons(bi,1)<y1 then, if y1<buttons(bi,2) then,
  705.          if buttons(bi,3)<y2 then, if y2<buttons(bi,4) then,
  706.             select bnames(bi)
  707.             case 'help' then,
  708.                helpme(1);
  709.                [i_i,y1,y2]=xclick();
  710.                hflag='on';
  711.             case 'stop' then,
  712.                write(%io(2),'Nothing to stop'),
  713.                [i_i,y1,y2]=xclick();
  714.                hflag='on';
  715.             case 'grill' then,
  716.                if gopt='on' then, 
  717.                   gopt='off';
  718.                else,
  719.                   gopt='on';
  720.                end,
  721.                makegrill(nx,nz,gopt);
  722.                [i_i,y1,y2]=xclick();
  723.                hflag='on';
  724.             case 'undo' then,
  725.                write(%io(2),'Nothing to undo'),
  726.                [i_i,y1,y2]=xclick();
  727.                hflag='on';
  728.             case 'quit' then,
  729.                qt=resume('on');
  730.             end,
  731.          end,end,
  732.          end,end,
  733.       end,
  734.    end,
  735.  
  736. //if line segment defined by first two points intersects a
  737. //previously defined line and y is in the frame then start a new line
  738.  
  739.       inflag='off';
  740.       if 1<=y1 then, if y1<=nx then,//y in the frame
  741.       if 1<=y2 then, if y2<=nz then,
  742.          inflag='on';
  743.          xsec=[];//list of intersections
  744.          nol=size(linelist);
  745.          for ln=1:nol,//test against all the lines for intersections
  746.             [testflag,bav]=testpt([x1;x2],[y1;y2],linelist(ln));
  747.             [br,bc]=size(bav);
  748.             xsec=[xsec,[bav;ln*ones(1,bc)]];
  749.          end,
  750.  
  751. //choose longest intersection which yields shortest segment
  752. //and add on the necessary extensions
  753.  
  754.          [xr,xc]=size(xsec);
  755.          if xc>0 then,
  756.             testflag='on';
  757.             btemp=xsec(1:2,:)-[x1;x2]*ones(1,xc);
  758.             [val,in]=maxi([1 1]*(btemp.*btemp));//longest intersection
  759.              xset("alufunction",6);
  760.              plot2d(x1',x2',[3,1],"000");//undo start circle
  761.              xset("alufunction",3);
  762.             x=xsec(1:2,in);//new x
  763.             x1=x(1);x2=x(2);
  764.             ln=xsec(4,in);li=xsec(3,in);
  765.             shortline=linelist(ln);//extension line
  766.             [sr,sc]=size(shortline);
  767.             if sc-li>li then,
  768.                yext=shortline(:,1:li);
  769.                [yer,yecf]=size(yext);//yec is useful for undoing
  770.             else,
  771.                yext=shortline(:,sc:-1:li+1);
  772.                [yer,yecf]=size(yext);
  773.             end,
  774.          end,
  775.       end,end,
  776.       end,end,
  777.       if inflag='off' then,//case where y is not in the frame
  778.              xset("alufunction",6);
  779.              plot2d(x1',x2',[3,1],"000");//undo start circle
  780.              xset("alufunction",3);
  781.          x1=y1;x2=y2;
  782.              plot2d(x1',x2',[3,1],"000");//undo start circle
  783.       end,
  784.       if testflag='off' then,//case where x-y does not cross a line
  785.              xset("alufunction",6);
  786.              plot2d(x1',x2',[3,1],"000");//undo start circle
  787.              xset("alufunction",3);
  788.          x1=y1;x2=y2;
  789.              plot2d(x1',x2',[3,1],"000");//undo start circle
  790.       end,
  791.    end,
  792.  
  793.    if flag='on' then,
  794.       plot2d([x1;y1],[x2;y2],[-1,1],"000"),
  795.       line=[yext,[x1;x2],[y1;y2]];
  796.       x1=y1;x2=y2;
  797.    end,
  798.  
  799. //continue line
  800.  
  801.    while flag='on',
  802.       flag='off';
  803.       [i_i,y1,y2]=xclick();
  804.  
  805. //check for a button
  806.    
  807.    [br,bc]=size(buttons);
  808.    hflag='on';
  809.    while hflag='on'
  810.       hflag='off';
  811.       for bi=1:br,
  812.          if buttons(bi,1)<y1 then, if y1<buttons(bi,2) then,
  813.          if buttons(bi,3)<y2 then, if y2<buttons(bi,4) then,
  814.             select bnames(bi)
  815.             case 'help' then,//get help
  816.                helpme(2);
  817.                [i_i,y1,y2]=xclick();
  818.                hflag='on';
  819.             case 'stop' then,//totally undo current line
  820.                [lr,lc]=size(line);
  821.                xset("alufunction",6);
  822.                plot2d(line(1,yecf+1:lc)',line(2,yecf+1:lc)',[-1,1],"000"),
  823.                xset("alufunction",3);
  824.                flag='off';
  825.                yecl=0;
  826.                line=0;
  827.             case 'grill' then,//toggle grill
  828.                if gopt='on' then, 
  829.                   gopt='off';
  830.                else,
  831.                   gopt='on';
  832.                end,
  833.                makegrill(nx,nz,gopt);
  834.                [i_i,y1,y2]=xclick();
  835.                hflag='on';
  836.             case 'undo' then,//undo line segment
  837.                [lr,lc]=size(line);
  838.                xset("alufunction",6);
  839.                plot2d(line(1,lc-1:lc)',line(2,lc-1:lc)',[-1,1],"000"),
  840.                xset("alufunction",3);
  841.                flag='off';
  842.                if lc>yecf+2 then,
  843.                   hflag='on';
  844.                   line=line(:,1:lc-1);
  845.                   x=line(:,lc-1);
  846.                   x1=x(1);x2=x(2);
  847.                   [i_i,y1,y2]=xclick();
  848.                   flag='on';
  849.                end,
  850.                if flag='off' then, yecl=0; line=0; end,
  851.              case 'quit' then
  852.                qt=resume('on')
  853.             end,
  854.          end,end,
  855.          end,end,
  856.       end,
  857.    end,
  858.  
  859.       if line<>0 then,//line drawing has not been stopped
  860.  
  861. //test if line intersects itself
  862.  
  863.          [lr,lc]=size(line);
  864.          if lc-1>yecf+1 then,
  865.             [testflag,bav]=testpt([x1;x2],[y1;y2],line(:,yecf+1:lc-1));
  866.          else,
  867.             testflag='off';
  868.          end,
  869.          if testflag='on' then,
  870.             write(%io(2),' '),
  871.             write(%io(2),'*********ERROR*********')
  872.             write(%io(2),' Lines are not allowed'),
  873.             write(%io(2),'to intersect themselves'),
  874.             write(%io(2),' '),
  875.             write(%io(2),'  Choose another point'),
  876.             flag='on';
  877.          else,
  878.  
  879. //test if point intersects a previously drawn line
  880.  
  881.          yext=[];
  882.          flag='on';
  883.          xsec=[];
  884.          nol=size(linelist);
  885.          for ln=1:nol,
  886.             [testflag,bav]=testpt([x1;x2],[y1;y2],linelist(ln));
  887.             [br,bc]=size(bav);
  888.             xsec=[xsec,[bav;ln*ones(1,bc)]];
  889.          end,
  890.  
  891. //check all the intersections and choose the shortest
  892.  
  893.          [br,bc]=size(xsec);
  894.          if bc>0 then,
  895.             flag='off';
  896.             btemp=xsec(1:2,:)-[x1;x2]*ones(1,bc);
  897.             [val,in]=mini([1 1]*(btemp.*btemp));
  898.             y1=xsec(1,in);y2=xsec(2,in);
  899.  
  900. //find the extension of the intersected line by the 
  901. //shortest segment of the intersecting line
  902.  
  903.             ln=xsec(4,in);li=xsec(3,in);
  904.             shortline=linelist(ln);
  905.             [sr,sc]=size(shortline);
  906.             if sc-li>li then,
  907.                yext=shortline(:,li:-1:1);
  908.                [yer,yecl]=size(yext);
  909.             else,
  910.                yext=shortline(:,li+1:sc);
  911.                [yer,yecl]=size(yext);
  912.             end,
  913.          end,
  914.  
  915. //plot line segment
  916.  
  917.          plot2d([x1;y1],[x2;y2],[-1,1],"000"),
  918.          line=[line,[y1;y2],yext];
  919.          x1=y1;x2=y2;
  920.       end,
  921.       end,
  922.    end,
  923.    if yecl<>0 then,
  924.       yec=[yec,[yecf;yecl]];
  925.    end,
  926.    write(%io(2),'Done With Current Line'),
  927.    
  928. //end
  929.  
  930. //[flag,bav]=testpt(p1,p2,line)
  931. //[flag,bav]=testpt(p1,p2,line)
  932. //macro which tests whether the line segment defined by the
  933. //two points p1 and p2 intersects any of the line segments
  934. //in line.
  935. // p1   :2x1 matrix giving the indices of point number one
  936. // p2   :2x1 matrix giving the indices of point number two
  937. // line :2xN matrix giving the indices of a line
  938. // flag :flag='on' indicates an intersection, flag='off' 
  939. //      :indicates no intersection
  940. // bav  :3xM matrix giving all the intersections found (if any)
  941. //      :and the position in the line of the intersection
  942. //
  943. //!
  944. // author: C. Bunks     date: 12-NOV-90
  945.  
  946. //set up arguments of fortran subprogram m45.f
  947.  
  948.    noi=0;
  949.    nc=maxi(size(line));
  950.    bav=0*ones(3,nc);
  951.    flag=0;
  952.    [flag,bav,noi]=fort('testpt',...
  953.                                    p1,1,'r',...
  954.                                    p2,2,'r',...
  955.                                    nc,3,'i',...
  956.                                    line,4,'r',...
  957.                                    flag,5,'i',...
  958.                                    noi,6,'i',...
  959.                                    bav,7,'r',...
  960.                                    'sort',...
  961.                                    [1,1],5,'i',...
  962.                                    [3,nc],7,'r',...
  963.                                    [1,1],6,'i');
  964.  
  965.    bav=bav(:,1:noi);
  966.    if flag=1 then, flag='on'; else, flag='off'; end,
  967. //end
  968.  
  969. //[]=redraw(linelist,seedlist,velolist)
  970. //[]=redraw(linelist[,seedlist[,velolist]])
  971. //Macro which redraws the lines drawn by velpic
  972. //and gives the locations and velocities chosen
  973. // linelist :list containing the lines to be drawn
  974. // seedlist :list containing the seed positions
  975. // velolist :list containing the velocity values
  976. //
  977. //!
  978. // author: C. Bunks     date: 12-NOV-90
  979.  
  980.    [lhs,rhs]=argn(0);
  981.  
  982. //determine the number of lines
  983.  
  984.    nol=size(linelist);
  985.    l1=linelist(1);
  986.    nr=maxi(l1(1,:));
  987.    nc=maxi(l1(2,:));
  988.  
  989. //draw frame
  990.     plot2d(0,0,[-1],"012",' ',[1,1,nr,nc]);
  991.  
  992. //draw the lines (the first line from velpic is the 
  993. //exterior line)
  994.  
  995.    for k=1:nol,
  996.       lk=linelist(k);
  997.       plot2d(lk(1,:)',lk(2,:)',[-1],"000"),
  998.    end,
  999.    if rhs=2 then,
  1000.       plot2d(seedlist(1,:)',seedlist(2,:)',[3,0],"000"),
  1001.    end,
  1002.    if rhs=3 then,
  1003.       [vr,vc]=size(velolist);
  1004.       toff=.05*maxi([nr,nc])/3;
  1005.       for k=1:vc,
  1006.          text=string(velolist(k));
  1007.          s1=seedlist(1,k);s2=seedlist(2,k);
  1008.          dess(50)=s1+toff;dess(52)=s2+toff;//positions of text
  1009.          plot2d(s1',s2',[3,0],"000");
  1010.          xstring(s1+toff,s2+toff,text,0,0);
  1011.       end,
  1012.    end,
  1013. //end
  1014.  
  1015. //[ind,indexlist]=id_rgn(indexlist,linelist,seed);
  1016. //[ind,indexlist]=id_region(indexlist,linelist,seed);
  1017. //Macro which determines all the indices of the matrix of dimension
  1018. //(nz X nx) which are in the region defined by the seed and the 
  1019. //linelist.  The elements which are in the same region as the seed
  1020. //are those which are on the same side of all the lines in the
  1021. //linelist.
  1022. // indexlist :(2xN) vector containing all the indices to be searched
  1023. // linelist  :list of lines
  1024. // seed      :pair of indices of the matrix identifying the region
  1025. // ind       :all indices of the matrix associated to the region
  1026. //           :defined by seed
  1027. //
  1028. //!
  1029. // author: C. Bunks     date: 12-NOV-90
  1030.  
  1031.    nlist=0*indexlist;
  1032.    ic=maxi(size(indexlist));
  1033.    nol=size(linelist);
  1034.    llist=[];noe=[1];
  1035.    for k=1:nol,
  1036.       noe=[noe,noe(k)+maxi(size(linelist(k)))];
  1037.       llist=[llist,linelist(k)];
  1038.    end,
  1039.    nolt=maxi(size(llist));
  1040.    noi=0;non=0;bav=0*ones(3,nolt);
  1041.    lln=0*ones(2,nolt);
  1042.    ind=0*indexlist;
  1043.    [ind,nlist,noi,non]=fort('id_rgn',...
  1044.                               indexlist,1,'i',...
  1045.                               llist,2,'r',...
  1046.                               seed,3,'r',...
  1047.                               ind,4,'i',...
  1048.                               nol,5,'i',...
  1049.                               nolt,6,'i',...
  1050.                               ic,7,'i',...
  1051.                               noi,8,'i',...
  1052.                               noe,9,'i',...
  1053.                               nlist,10,'i',...
  1054.                               bav,11,'r',...
  1055.                               lln,12,'r',...
  1056.                               non,13,'i',...
  1057.                               'sort',...
  1058.                               [2,ic],4,'i',...
  1059.                               [2,ic],10,'i',...
  1060.                               [1,1],8,'i',...
  1061.                               [1,1],13,'i');
  1062.    ind=ind(:,1:noi);
  1063.    indexlist=nlist(:,1:non);
  1064.  
  1065. //end
  1066.  
  1067.  
  1068. //[]=makegrill(nx,nz,gopt)
  1069. //Plot a grill in the region [1 nx] x [1 nz]
  1070. // nx   :Number of x positions
  1071. // nz   :Number of z positions
  1072. // gopt :Plotting flag where 'on' means plot 
  1073. //      :and 'off means unplot
  1074. //
  1075. //!
  1076. // author: C. Bunks     date: 12-NOV-90
  1077. // Change JPC 2 mars 1992
  1078. //   dess;
  1079.    if gopt <>'on' then,xset("alufunction",6);end
  1080.    for k=2:nx-1, plot2d([k;k],[1;nz],[-2],"000"), end,
  1081.    for k=2:nz-1, plot2d([1;nx],[k;k],[-2],"000"), end
  1082.    xset("alufunction",3);
  1083. //end
  1084.  
  1085. //[vi]=velcalc(indexlist,seedlist,velolist)
  1086. //Create velocity field for a region defined by indexlist
  1087. //where velocity varies linearly between control points
  1088. //defined by the seedlist sl and its associated velocities
  1089. //in vl.
  1090. // indexlist :list of indices of region
  1091. // seedlist  :indices of velocity seed locations
  1092. // velolist  :velocity control points
  1093. // vi        :velocities for region
  1094. //
  1095. //!
  1096. // author: C. Bunks     date: 12-NOV-90
  1097.  
  1098.    [vr,vc]=size(velolist);
  1099.    [ir,ic]=size(indexlist);
  1100.    vi=velolist(1)*ones(1,ic);
  1101.    for k=1:vc-1,
  1102.       s1=seedlist(1:2,k);s2=seedlist(1:2,k+1);
  1103.       v1=velolist(k);v2=velolist(k+1);
  1104.       x=indexlist(1,:);y=indexlist(2,:);
  1105.       x1=s1(1);y1=s1(2);
  1106.       x2=s2(1);y2=s2(2);
  1107.       dx=x2-x1;dy=y2-y1;dv=v2-v1;
  1108.       r2=dx*dx+dy*dy;
  1109.       xv=x-x1*ones(x);
  1110.       yv=y-y1*ones(y);
  1111.       gr=(xv*dx+yv*dy)/r2;
  1112.       vi=vi+mini(maxi(gr,0*gr),ones(gr))*dv;
  1113.    end,
  1114.  
  1115. //end
  1116.